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Oscillations are ubiquitous phenomena in biological systems. Conventional models of biological 
periodic oscillations usually invoke interconnecting transcriptional feedback loops. Some specific 
proteins function as transcription factors, which in turn negatively regulate the expression of the 
genes that encode these "clock proteins" . These loops may lead to rhythmic changes in gene expres- 
sion in a cell. In the case of multi-cellular tissue, collective oscillation is often due to synchronization 
of these cells, which manifest themselves as autonomous oscillators. In contrast, we propose here 
a different scenario for the occurrence of collective oscillation in a group of non- oscillatory cells. 
Neither periodic external stimulation nor pacemaker cells with intrinsically oscillator are included 
in present system. By adopting a spatially inhomogeneous active factor, we observe and analyze a 
coupling- induced oscillation, inherent to the phenomenon of wave propagation due to intracellular 
communication. 

PACS numbers: A PACS will be appear here 
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I. INTRODUCTION 

Oscillation is ubiquitous in nature, not only in physics 
and chemistry but also biology. Biological oscilla- 
tions can be observed over a wide range of time- and 
population-scales, from a circadian rhythm of about 24 
hours [H to a segmentation clock of less than 2 hours , 
from whole-body oscillatory fevers Q to periodic protein 
production in a single cell [4]. On the other hand, sound 
theoretical studies have been undergoing since long be- 
fore the observation became possible in molecular level. 
There are many theoretical models to explain these phe- 
nomena. Despite their diversity of biological insights, 
these models share some common points. 

Proteins are produced by the transcription and trans- 
lation of specific sequences of DNA. On the other hand, 
proteins can bind to a transcription promotor on DNA 
and hence suppress or enhance gene expression. A tran- 
scriptional negative feedback loop [1, Q and a delay 0, [1| 
in the inner cellular gene-protein network are considered 
to be important elements that contribute to the oscilla- 
tory expression of DNA and protein production. From 
the perspectives of dynamical systems, such oscillations 
are limit cycles that can be generated from Hopf bifur- 
cation by choosing an appropriate parameter set and ini- 
tial condition. Consequently, in the case of a cell group 
or multi-cellular organism with an oscillatory character, 
such as cardiac tissue and a segmentation clock in the tail 
of PSM (Presomitic Mesoderm), the synchronization of 
coupled oscillators is often used to explain the observed 
collective oscillation [l^, |Tl[ . 

However, periodic oscillation is only a small part of 
the dynamical behavior of a cell. Oscillation may cease 
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if the conditions are changed, and most cells tend to settle 
into a seeming stable state. For example, electrical ac- 
tivity in /3-cells exhibits slow periodic oscillation at the 
macro scale of islets of Langerhans, while much faster 
excitability instead of oscillation when isolated [H, ■ 
In another example, the three proteins (KaiA, KaiB, and 
KaiC) , identified as important for the circadian rhythms 
in cyanobacterium Synechococcus elongates, behave as a 
bistable toggle switch due to a double-negative-feedback 
loop. Oscillation could then arise from the successive 
switch between these two stable steady-states [13, [T5| . 
Moreover, most recent studies also suggested that nega- 
tive transcriptional feedback is not sufficient, and in some 
cases not even necessary, for circadian oscillation. In- 
stead, intracellular signaling, such as that involving Ca^+ 
and cAMP, together with transcriptional feedback plays 
a key role in long-term circadian pacemaking . These 
evidences raise the possibility that intrinsic oscillatory 
cells are not indispensable in an oscillatory organism. 



In this paper, we study the occurrence of collective 
oscillation from non-oscillatory system. In contrast to 
the conventional mechanism of synchronized oscillators, 
none of the individual cells in our model is intrinsically 
oscillatory. A few studies in the context of mathemat- 
ics and physics have revealed the possibility of collec- 
tive oscillating patterns. The first example was pro- 
posed by Smale [13], who found that two "dead" cells 
can become "alive" via diffusive coupling. More re- 
cently, other studies have examined this behavior in de- 
tail [3i [I3l ■ In-phase and anti-phase self-sustained os- 
cillation of excitable membrane via bulk coupling have 
been observed [2^1 ■ The models considered in these re- 
ports have mostly involved coupled identical excitable 
cells with mono-stability. Some more complicated ap- 
proaches include, for example, using a unidirectional cou- 
pling scheme [21,], applying a periodic stimulation 
coupling the system with an oscillatory boundary 
introducing heterogeneity into excitable media [2J 
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tivity propagating in discrete cellular automata model 
[25| . and so forth. A commonly used idea is to set iso- 
lated cells at a subthreshold quiescent state, and then 
push them over into the oscillatory regime to generate 
pacemaker cells by extra force or coupling. That means 
cells are possible to manifest themselves as oscillators. 
However, little attention has been paid to the emergence 
of oscillation in systems that are completely independent 
of oscillating elements. Unlike previous studies, geomet- 
rical structure of nuUclines of cells in our model prevent 
dynamics from being oscillation. There are two "engines" 
in our model to drive the self-sustained collective oscilla- 
tion, neither of which is oscillatory pacemaker. The one 
is bistable cell switching between two stable states, the 
other is mono-stable cell with excitability. Two engines 
work cooperatively due to the wave propagation. 

On the other hand, in a bistable system, a stationary 
front can bifurcate into a pair of fronts that propagate in 
opposite directions, which is known as non-equilibrium 
Ising-Bloch (NIB) bifurcation [2i,[27|. Perturbation for 
the occurrence of NIB bifurcation can be induced by lo- 
cal spatial inhomogeneity [1^. A more global analysis 
showed that the NIB point is only part of the story, and 
concluded that an unstable wave front is intrinsic to me- 
dia that are spatially inhomogeneous [H, [sOj- An un- 
stable wave front may manifest itself as a reflected front, 
tango wave [3l|, pacemaker |3^] and so on. In this paper, 
we think about these phenomena beyond mathematics 
and physics, and extend their application to biological 
oscillators. 

Moreover, although most studies have been performed 
on a spatial continuum described by partial differential 
equations (PDE), continuum models neglect the effects 
of cellular discreteness 33]. In fact, from the viewpoint 
of biology, the size of cells can not decreased infinitely. 
This intrinsic property is difficult to ignore, especially at 
the stage of initial development of an organism, when the 
cell size is comparable to that of tissue. In addition, there 
are mathematical reasons to explore the system dynamics 
with spatial discretization. PDE and ODE (ordinary dif- 
ferential equations) have different theoretical frameworks 
and produce different results. Several significant features 
of discreteness, such as wave propagation failure [s^ . can 
not occur in a continuum model. Therefore, in this paper 
we will consider an array of spatially discrete cells, and 
discuss the impact of discreteness. 



II. DESCRIPTION OF THE MODEL 

A. One-dimensional cellular array 

In this paper, we consider cells in one-dimensional 
space. Cells are coupled by intracellular signaling 
molecules, which flow through channels in a mem- 
brane due to concentration difference or depolarization- 
mediated flux. The intracellular signaling small-molecule 
can be produced by a series of process from some genes 



functioned as activator, and then trigger transcriptional 
feedback loops of adjacent cells. We assume that the cou- 
pling interaction takes place in a diffusion-like manner. 
If we include an inhibitor, which can locally repress the 
expression of activator genes, a one-dimensional array of 
N cells can be described as 



Ui = f{ui,Vi, Pi) + D{u^_i + u^+i 



2m0 



(1) 
(2) 



where u and v are concentration of activator and in- 
hibitor, respectively, i G {l...A^} is the index of the 
cell in the chain and D is the coupling strength of u. / 
and g are corresponding reaction functions. The bound- 
ary condition is zero flux, i.e., uq = ui and un = un+i- 
Finally, Pi is an environmental parameter, which will be 
discussed in detail later. 

In this study, we only consider coupling of the activa- 
tor. For most of the models that have been used to study 
pattern formation, diffusion is assumed to occur for ev- 
ery elements. Specially, much greater diffusion of the 
inhibitor is necessary to induce Turing instability [ssj . 
However, a cells membrane is very selective for passage 
of substances. Complicated intracellular reactions usu- 
ally take place locally, but are triggered by only one or 
a few specific signaling molecules. For example, while 
the segmentation clock involves the cyclic expression of 
many genes, the crucial pathway for coupling only in- 
volves the transmembrane receptor Notchl [36|. Thus, 
in the context of biology, we only consider coupling with 
the activator, and the inhibitor in our model is merely a 
local state variable. 



B. Active factor 

The development of a multicellular organism begins 
with a single cell, which divides and gives rise to cells 
with different typologies. Different cells are organized ac- 
cording to certain secreted chemicals, called morphogens. 
Despite improvements in experimental and theoretical 
approaches, the mechanisms of morphogenesis are still 
unclear. Usually, morphogens are considered to be pro- 
duced at specific sites and diffuse through the organ- 
ism [37[ . Quite recently, evidence of a "shuttling-based" 
mechanism has been presented f38j. The key in such 



models is their ability to define a robust and scaling 
profile, usually a concentration gradient, of morphogens. 
More broadly, we can suppose that some environmen- 
tal parameters act as morphogens. The environment in 
which an organism develops supplies nutrition for growth, 
and the intracellular volume in direct contact with the 
border gets more and that deep inside cells gets less. 

In this paper, we do not consider any specific chemi- 
cal substance, and instead merely suppose that there is 
a certain factor, which we refer to as the active factor, 
to obtain information regarding the relation between po- 
sition and cell dynamics. The above-mentioned active 



FIG. 1: Profile of Ti obtained from Eq.©, wlien To = 15, ^ = 
5. 

factors can affect the fate of cells in a concentration de- 
pendent manner [39l |40| . 

Without losing generality, we assume that the active 
factor r is constant at the boundaries of an organism, 
where the source site of morphogens are usually located. 
It diffuses into the organism field with a diffusion con- 
stant Da, and is degraded at rate a. Thus, we have 

a^-^«a^-"r. (3) 

Since our model is based on coupled ODEs independent 
of spatial variation, the profile of the active factor satis- 
fies a scaling property. By normalizing the field size to 
one, we can get a steady profile {dV/dt = 0) of F as 

r(x) = _iiL^((e-« - l)e«- - (e« - l)e-«-), (4) 

where Fq = F(0) = F(l) is the value at two boundaries, 
and f = 1/A = ^Ja/Da is the inverse of the decay length. 
A typical profile of T{x) is shown in Fig. [TJ Circles in- 
dicate the value of F for discrete cells {N = 20 in the 
figure) placed uniformly in the scaling field. 

C. Bistability 

We assume that cells normally prefer to live in a stable 
state, and cells with a high concentration of active factor 
are capable of switching between two states. This kind 
of bistability is very important and has been observed 
in various biological systems [4l[. For example, the ex- 
pression of the Dictyostelium cAMP phosphodiesterase 
gene behave as a bistable switch employing intracellular 
cAMP as a regulator of cell fate [4^ , the Cdc2 activation 
system in Xenopus egg extracts is bistable and charac- 
terized by biochemical hysteresis |43| . the inducible lac 
operon in E. coli shows bistability '44'| , and so on. 

Usually, bistability arises from positivejOr double- 
negative genetic regulation loops [4l|, |4^ |4^. It was 
recently suggested that stochastic fluctuation plays an 
important role in the nature of the transition between 
bistable states [U Ell. Moreover, physical regulation 
of protein production, which has been much less con- 
sidered by biochemists, also plays an important role in 
the origin of the bistability. It has been observed that 
discrete transition between folding and unfolding states. 



(a) (b) 

FIG. 2: Nullcline diagrams in bistability (a) and mono- 
stability (b), respectively. 

namely a first-order phase transition, can take place in 
giant DNA [49[. Similar discrete switch can also occur 
for RNA ^5^, protein [5i| and other molecules 52). This 
discrete transition leads to the ON/OFF switching of the 
production of a specific protein. 

D. Model equations 

We describe the dynamical reaction function of each 
cell by using the two-component Fitzhugh-Nagumo equa- 
tions 

u = fiu, V, F) = Tu(u — a){l — u) — V (5) 
V = g{u, v) — e(/3w — v) (6) 

where m is a variable related to the expression level of 
specific activator genes, v is the inhibitor to repress u, 
e, which is much smaller than 1, is the slower growth 
factor of inhibitor v and F is the active factor discussed 
previously. Note that the kinetics of inhibitor here is a 
rather natural unit process in many of biochemical reac- 
tions. Throughout this paper, the following parameters 
are fixed 

a = 0.3, /3 = 0.5, e = 0.02. (7) 

The Fitzhugh-Nagumo model has been well studied for 
description of excitable behavior in biology. Rich non- 
linear dynamics can be observed by tuning parameters. 
Specifically, with the above parameters, the model is 
mono-stable at small value of F, and will happen a saddle- 
node bifurcation at F = 4.08 and a Hopf bifurcation at 
F — 4.27, which leads to bistability. Thus, in the case 
of the spatial profile of F^ as shown in Fig. [H only the 8 
central cells (7*'* - 14*'*) are mono-stable, while the others 
are bistable {Tq = F15 = 4.37). 

Figures [2la)&(b) show the nuUclines with bistability 
and mono-stability for when F is large and small, respec- 
tively. Again, none of the cells show oscillation in the 
absence of coupling. More important, from the geometry 
property of nuUclines in the figure, no oscillatory condi- 
tion could be found by moving cubic nullcline (u — 0) up 
and down. That is, no pacemaker cells can be generated 
from activator coupling. 
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FIG. 3: Collective oscillation observed in a chain of cells when 
D = 0.7 in Eqs. ([8l?i [H. (a) Spatio-temporal plot of the 
collective oscillation of Ui . The black and white indicate Ui — 
1 and Ui = 0, respectively, (b, c) Waveforms of u (solid) and 
V (dashed) in the 3rd and 4th cells. 

If we substitute Eq. dSMH) into Eq. dHyU, we get the 
system equations used in this paper. 

iii = riUi{u, - a){l ~ u^) - Vi (8) 

+ D{ui-i + Ui+i - 2ui) 
Vi = e{Put - Vi) (9) 

When we change the couphng strength D, we observe the 
occurrence, variation and disappearance of self-sustained 
collective oscillation in the cell array. 



III. SELF-SUSTAINED COLLECTIVE 
OSCILLATION 

A. Normal collective oscillation 

Figure ^a.) shows a typical oscillation when the cou- 
pling strength D — 0.7. Figure [IJa) shows an enlarged 
view of a single period of oscillation. As the initial condi- 
tion, we set the 1st cell as being excited, since stimulation 



is usually input from the border. Initially, (0 < t < 80), 
a traveling wave appears due to excitation at the border. 
The traveling front then sweeps over the cell array and 
makes all of the cells excited (see Fig. UJ^b)). Although 
the central cells are also turned ON due to the interaction 
with other cells, they can not stay in the excitable state 
for a long time. Instead, they soon return to their sta- 
ble equilibrium (see Fig. UJ^c)), and hence generate two 
counterpropagating wave backs, as shown in Fig. HJ^d). 
These two wave backs propagate outward until the 3'"'^ 
and 18"* cells and stop suddenly due to the spatial dis- 
creteness (see Fig. HJe)). The "wall" cells do not jump 
from the ON state to the OFF state and only exhibit 
slight oscillation closed to their equilibrium. As an exam- 
ple, the difference between the 3''^ and 4*'* cells is shown 
in Fig. ^h, c). At this critical interface, the inhibitor 
V slowly decreases so that the 4*'* and 17*'' cells restore 
excitability after a while. The central cells can then be 
excited again by the pair of reflecting wave fronts, as 
shown in Fig. [Iff). Pushed by the wave, the central cells 
will be excited again. This process repeats and causes 
the collective oscillation inside multi-cell tissue without 
oscillatory cells. 



B. Stationary state before birth of oscillation 

The above collective oscillation can be observed when 
the coupling strength is larger than a threshold, below 
which wave backs (see Fig. IDJd)) fail to reflect, and the 
state in Fig. UJe) is maintained. Figure [5] shows a spatio- 
temporal diagram, where the central cells stay silent 
while excited bands appear close to the two borders. 

Note that this phenomenon could not take place in 
a continuum counterpart. The existence of a coupling 
strength threshold under which wave propagation failure 
occurs is unique to a spatially discrete system. In addi- 
tion, there is another threshold, which is even smaller, 
for which the wave front stops propagating. In this case, 
the excited signal at the border fails to propagate for- 
ward, but we would like to postpone this interesting phe- 
nomenon on another paper, since it is less related to 
present work. 



C. In-phase and anti-phase period doubling 
oscillation closed to the boundary 

With an increase in the coupling strength D, the char- 
acteristics of oscillation can be changed. Figure [6] shows 
that the position of oscillation periodically shifts. The 
3'"'' and 18"" cells oscillate with a nearly doubled period, 
in anti-phase (Fig. [HIb,c)). Globally, tissue oscillates in 
two groups with the same cell populations but different 
positions: No. 3-No. 17 (15 cells) and No. 4-No. 18 (15 
cells), respectively. 

Interestingly, by slightly increasing the coupling 
strength D, say to D = 0.9, we found a different type 
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FIG. 4: (a) Spatio-temporal diagram of Ui over a single pe- 
riod, (b)-(f) Snapshots of u and v at several time points in 
one period, where the horizontal axis is cell number from 1 to 
20. This illustrates the change in wave propagation at differ- 
ent stages. Solid curves and dashed curves indicate u and v, 
respectively. An animation, through which the behavior can 
be understood more intuitively, is available at [s^ ]. 



of period doubling, as shown in Fig. [T] For comparison 
with the case oi D = 0.8, although the critical interface 
between ON and OFF shifts periodically as in Fig. [6l 
there is no phase difference between the 3'"'' and 18*'' 
cells. As is clearly shown in their waveform (Fig.[7I^b,c)), 
these two boundary cells oscillate in-phase, instead of 
anti-phase (Fig. [6fb,c)). Therefore, in the present condi- 
tion, a periodic change does not take place in the position 
of oscillation. Instead, the population of oscillating cells 
changes. More precisely, tissue oscillates in two groups: 
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FIG. 5: Spatio-temporal diagram of Ui in a stationary state. 
Wave propagation stops and no oscillation occurs in the case 
of weak coupling D = 0.47. 
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FIG. 6: Anti-phase mode in period doubling produces collec- 
tive oscillation with a periodic position shift when D — 0.8. 
(a) Spatio-temporal diagram of Ui. The grayscale black and 
white indicate Ui — 1 and Ui = 0, respectively, (b) and (c) are 
waveform diagrams of the 3rd and 18th cells. Activator u and 
inhibitor v are shown in solid and dashed curves, respectively. 



No. 3-No. 18 (16 cells) and No. 4-No. 17 (14 cells), re- 
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FIG. 7: In-phase mode in period doubling produces collec- 
tive oscillation with a change in the periodic population when 
D — 0.9. (a) Spatio-temporal diagram of Ui. The grayscale 
black and white indicate Mi = 1 and Ui = 0, respectively, (b) 
and (c) are waveform diagrams of the 3'''* and IS"' cells. Acti- 
vator u and inhibitor v are shown in solid and dashed curves, 
respectively. 

spectively. 

Moreover, by setting the initial condition of the cells 
identically, i.e., all in the ON state at i = 0, we found 
checked that same symmetric collective oscillation can 
also occur in the case oi D = 0.8. Therefore, we con- 
clude that these two types of oscillation are caused by 
the same bifurcation. Because the wavefront propagates 
faster with larger D, a larger coupling strength can re- 
duce the time lag between the two boundary cells being 
stimulated. If the time lag is smaller, the two boundaries 
converge to in-phase oscillation. On the other hand, if 
the time lag is large, they will exhibit anti-phase oscilla- 
tion. 



D. Oscillation death 

With a increase in coupling strength D, we observed 
that the change in the periodic position or population 
stopped, and normal oscillation returned. In comparison 
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FIG. 8: Spatio-temporal diagram of Ui. Black and white 
indicate m — 1 and Ui — 0, respectively, (a) D — 2.53, 
oscillation starts to collapse; (b) D — 2.6, oscillation ceases 
after one cycle. 



to the case of 13 = 0.7, the total population of oscillating 
cells increased from 14 (No. 4 to No. 17) to 16 (No. 3 to 
No. 18). 

The oscillation suddenly dies when D is as large as 2.6. 
Figure El^b) clearly shows that the central cells start to 
oscillate after all of the cells are excited, but this oscilla- 
tion is not sustained. In this strong coupling condition, 
the boundary cells can not recover their excitability, so 
that the wave front propagating from the center is unable 
to stop and reflect to generate successive oscillation. 

Before the oscillation stops, there is a narrow parame- 
ter region of 2.53 ^ D ^ 2.56, where only one side of the 
"wall" alternatively collapses, and a complicated period- 
4 collective oscillation is observed (Fig. [5Ka)). 

E. Overall perspective and bifurcation 

There are many factors that may influence the oscil- 
lation behavior. For example, if the spatial profile of 
active factor F becomes more "steep" rather than a gen- 
tle slope, the wave back will tend to be locked and fail 
to reflect from boundary. Moreover, if we reduce the 
excitability of cells by increasing e, it will be more dif- 
ficult for wave front to propagate cross the center, and 
only the half part with stimulation can oscillate. Since 
the global oscillation is induced by mutual coupling, we 
are going to study the oscillation behavior with respect 
to the coupling strength. Here, we sweep D from 0.45 to 
2.7, and summarize the variation in the oscillation period 
and position of the left border of oscillation region. 

If the coupling strength is smaller than 0.48, there is 
no oscillation, and 4 cells from the tissue border are ex- 
cited while cells 5-16 are silent. Oscillation takes place 
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FIG. 9: Phase diagram of (a) cell number for the left bound- 
ary of collective oscillation, (b) the oscillation period, (c) an 
enlarged view with in-phase and anti-phase period two solu- 
tions, with respect to the change in the strength of coupling. 

when the wave back passes the 4*'' cell aX D = 0.48. The 
border then shifts between 3 and 4, while in-phase and 
anti-phase period doubled oscillation occur, roughly be- 
tween 0.78 < D < 0.97. Finally, the oscillation reaches a 
maximum region: from cells 3 to 18, until D is too large 
for oscillation to occur. Figure [nja) shows the expansion 
of oscillatory region. 

Variations in the period of oscillation are shown in 
Fig-IHb). Once the central cells start to collectively os- 
cillate. The period rapidly decreases when the coupling 
strength increases. The rate of the period decrease grad- 
ually slows. The period changes little in the region where 
D is large. This phenomenon occurs because the station- 
ary interval (Fig.dfe)) greatly contributes to the period 
of oscillation. The decrease in the stationary interval 



significantly shortens the period of oscillation when D 
is small. However, when D is large enough, the wave 
backs reflect immediately without stopping, and the pe- 
riod is determined mainly by the velocity of propagation. 
Therefore, the presented oscillation is robust at strong 
coupling condition, and tunable at weak coupling case. 

There is a parameter region (the curve of period is 
drawn in dashed curve Fig. [9Kb)) in which system under- 
goes Period-Doubling bifurcation and the period- 1 solu- 
tion lose its stability. Meanwhile, period two solutions 
appear around this region. We show more details in 
Fig-EUc). In the figure, we draw two intervals in a period 
two solution, by measuring the time when M4 positively 
cross the section: U4 = 0.5. Circles (o) and crosses (x) 
indicate in-phase and anti- phase solutions, respectively. 
The in-phase period-2 solution is the result of Period- 
Doubing bifurcation, while the anti-phase period-2 so- 
lution occurs saddle-node bifurcation. The anti-phase 
period-2 solution has a wider parameter region than in- 
phase one, and coexistence with fundamental period-1 so- 
lution can be observed in both side of D. There are quite 
complicated bifurcation phenomena specially around the 
occurrence of Period-Doubling bifurcation. We have even 
found period-3 solution (in-phase one around Z)=0.802, 
and anti-phase one around D = 0.975, respectively). Al- 
though they are very interesting in the viewpoint of non- 
linear dynamical system, we leave them to our future 
work, because current paper is going to discuss the pos- 
sibility of global oscillation and its potential applications. 

IV. DISCUSSION 
A. Discreteness vs. continuum 

The above phenomena are observed in a spatially dis- 
crete system, described by ordinary differential equa- 
tions. As briefly introduced in Sec[Tl this discreteness is 
important in both a mathematical and biological sense. 
Let us discuss this significance in more detail. 

The diffusion term D{d^u/dx^) in a one-dimensional 
spatially continuous reaction-diffusion model can be for- 
mulated as D{ui^i + Ui-i-i — 2ui)/ IS.x^ in its difference 
version. This type of conversion is a common approach 
to solving PDE numerically. The diffusion rate D usually 
does not change much for a specific substance under con- 
stant conditions. Thus, if we assume that the coupling 
is mainly due to the diffusion-like effects of substances, 
the coupling strength D « I?/Aa;^ changes in square or- 
der with respect to variation of Aa:, which biologically 
corresponds to the distance between cells or the cell size. 
Since the profile of the active factor has a scaling prop- 
erty, it is reasonable to suppose that this gradient works 
for a field of any size. Thus, we can study how a change 
in number of cells N and distant of cells Ax affects global 
dynamics. 

Figure llOf a) shows spatio-temporal diagrams with a 
10-fold increase in the number of cells, N = 200. Other 
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FIG. 10: Spatio-temporal diagram of m. Black and white 
indicate Ui = 1 and Ui = 0, respectively. A'^ = 200. (a) 
D = D/Ax'^ = 0.7. (b) D = D/Ax'^ = 70. 



parameters are the same as those in Fig. [3l Obviously, 
more time is required for a wave to sweep over the organ- 
ism. The period of oscillation and the phase difference 
between the two sides increase greatly. On the other 
hand, if the distance Ax between cells becomes smaller 
and smaller when cell population increases, the system 
manifests itself more like a continuum than a discrete 
system. In this case, the coupling strength will increase 
dramatically as a square with respect to the decrease in 
Ax, i.e., a 100-fold increase in D in present case. When 
D is as large as 70, we have the spatio-temporal diagram 
given in Fig. fTOTb) . When we compare this with Fig. [3l 
there is little change in the period of collective oscillation. 
This suggests that the clock tends to run more punctu- 
ally. In a mathematical sense, when the population of 
cells is large enough in a fixed field, the behavior of the 
organism will follow the solution of a specific partial dif- 
ferential equation, which is independent of the number 
of cells. 

Figure. [TOT a) simply corresponds to the case that cells 
grow in an open space, and extend the field by keeping 
size and distant of cells constantly. During the initial 
period of development, however, cell divisions within the 
egg proceed quickly, without much increase of the total 
cell mass and size. Thus, cells at this stage rapidly de- 
crease in diameter. This may interpret biologically the 
situation of Fig. fTUfb). 

There are many biological situations, however, that the 
intercellular coupling does not follow the diffusion-like 
ways, such as communication involving the delta-notch 
signaling pathway [s^ . In those cases, Aa; has few direct 
infiuences on D, which may represent "bottlenecks" ir- 
respective to the diffusion. Thus, the modeling based a 
continuum is inappropriate for some conditions. 




FIG. 11: Phase portrait diagrams with snapshots of the dy- 
namical nuUcline. Rows indicate the time evolution from the 
top down, and columns indicate the number of cells (3 at left, 
4 at middle and 5 at right). Dashed curves are the limit cycle 
solution. Cubic function curves are the nuUcline of Mi — 0. 
Straight lines are the nuUcline of Vi — 0. Circles are the 
position of {ui,Vi) at specific times. An animation of the dy- 
namical nuUclines can be found at 1541. 



B. Understanding the mechanism 

Self-sustained collective oscillation is caused by the 
excitability of cells and their mutual interaction. The 
system involves complicated bifurcations. We present 
here some qualitative ideas regarding how this oscillation 
takes place. 

From dynamical equations (j5l6p and their nuUcline 
shown in Fig. [21 we know that a single cell can exhibit 
either bistability or mono-stability. However, if we in- 
troduce coupling, the geometry of nuUclines of one cell 
will dynamically change according to its own state and 
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those of its neighbors. Because it was assumed that the 
communication between cells is only mediated via the 
activator u, the strict nuUcline Vi = is independent of 
coupling. 

V, = G(uO = f3u,. (10) 

From Eqs.([TB45|), we obtain the function for the nuUcline 

lii = as 

Vi = F{ui,ri) = TtUi{ui - a)(l - Ml) + AUi, (11) 

where AUi = D{ui+i + ui^i — 2uj) is the offset of the 
cubic function due to coupling. Thus, the nuUcline v = 
F{-) dynamically moves up and down in the phase plane, 
corresponding to the state of Ui-i, Ui, Ui+i. 

In Fig. [TI] we show the phase portrait of cells around 
the oscillation border (cells No.3-5), as well as their dy- 
namical nuUcline at some turning points. Snapshots are 
taken under the same conditions of normal oscillation as 
shown in Fig. [3l 

The first row, (a, b, c) of Fig. [Tl] are all in the ex- 
cited state, i.e., for all three cells, Ui is close to 1. 
Therefore, under this condition, the vertical offset AC/ 
of the nuUcline is nearly zero, and all three cells exhibit 
bistability. The second row is taken at t = 226, when 
the wave back comes (see Fig. Iljd)), and the 5*'' cell 
moves towards its lower equilibrium (Fig. Illf f)). Since 
AUi — D{u3 -f Ms — 2M4), a sudden drop in M5 leads 
to a rapid decrease in the cubic nuUcline. As shown in 
Fig. fTlTe). the cubic nuUcline moves down so that the 
higher equilibrium disappears. Thus, the 4*^ cell be- 
comes mono-stable and the state quickly converges to 
the left branch of the cubic nuUcline. The decrease in M4 
pushes At/4 back to a positive value, and makes the nuU- 
cline of the 3'"'' cell move down, as shown in Fig. fTTT g. h). 
However, since the 3'"'' cell has a larger F, which controls 
the amplitude of the cubic nuUcline, even if M4 decreases 
to its lowest value (Fig. [TlT h)). i.e., AU3 reduces to its 
minimum, the cubic and straight nuUclines still intersect, 
and the higher equilibrium remain. This explains why the 
wave back passes the 4*'' cell, but stops at the 3'"'' cell 
(Fig. HKa)). After propagation stops, there is a relatively 
long refractory period from time 230 to 300. In this in- 
terval, there is a slow decrease in the inhibitor W4. Since 
M3 ^ 0.8 and M5 w 0, although slight increase in M4 pulls 
the cubic nuUcline down, AC/4 is still so large that the cu- 
bic nuUcline is above the straight nuUcline (Fig. [Tljk)). 
Under this condition, the cell is mono-stable, with the 
equilibrium at the right branch of the cubic function. 
Thus, after a while, the state of u will switch to a higher 
value (Fig. [TlT n)). and leads to a reflecting wave front 
(Fig. m^f)). Finally, the states return to the situation of 
Fig. Illf b) after another refractory period. 

Note that a smaller coupling strength D leads to a 
smaller offset AC/. If we move down the cubic nuUcline 
slightly to cross the straight nuUcline in Fig. [TTJk), the 
4*'^ cell becomes bistable. This will disable the switch 
from left to right, and stop the oscillation (Fig. (5]). 



From the above description and Fig. [Tl] we conclude 
that the boundary cell, here No. 4, which is bistable with- 
out coupling, turns to switch between two types of mono- 
stable dynamics. As introduced in Sec. HI it is different 
from the studies changing nuUclines via coupling to an 
oscillatory geometry. This switching becomes the power 
that underlies the self-sustained oscillation observed in 
the present model. The variation of the offset of the dy- 
namical nuUcline of the boundary cell gives rise to rich 
oscillation phenomena. 

C. Conditions for oscillation 

We will now explore the conditions for oscillation in 
an approximate manner by studying the dynamics on 
the oscillation border, where wave backs (WB) stop and 
wave fronts (WF) generate. Based on an investigation of 
the dynamical nuUcline and state variable, we concluded 
that a wave back will not pass a critical cell, c, if the 
nuUclines still intersect at the right branch when cell c+1 
has dropped to its lower equilibrium (Fig. Illf g.h)). In 
contrast, if the intersections disappear, the state of Uc 
will switch to a lower equilibrium. Then cell c is possible 
to oscillate, if it fires a wave front, in the case that the two 
nuUclines do not cross at the left branch when the state 
of the inhibitor recovers to its lower limit (Fig. [TIT e.f)). 

Thus, we can roughly solve the condition by finding 
two possible tangency for the two nuUclines Eq. ^TU\i and 
Eq. pT|) . This can be achieved using the following equa- 
tions: 

di^(M,r) _ dG(M) 



Equation is for two nuUclines with the same slope. 
By substituting F and G into Eq. (fT^ . we have 

T{-3u^ + 2{l + a)u-a)^P, (13) 

from which we obtain two solutions 

For cell c to propagate a wave back, there should 
be only a lower equilibrium when Uc close to the 
higher tangency point. The corresponding condition is 
F{uti,Tc) < G{uTi)- Substitution leads to 

rcUTi(MTi ~a){l-UTj + D{Uc-l + Uc+i-2uTi) < [iuT^- 

(15) 

On the other hand, for cell c to generate a wave front, 
there should be no lower equilibrium when Mc is close 
to the lower tangency point. This simply means that 
F{ut2,^ c) > G{ut2), which can be rewritten as 

rcUT2(uT2 - a){l - UT2) + D{Uc-l + Uc+l - 2UT2) > (3ut2- 

(16) 

Two critical conditions are shown in Fig. 1121 
Approximation 
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FIG. 12: Schematic diagram of two critical tangency situa- 
tions, corresponding to the conditions for which (a) a wave 
back passes and (b) a wave front is generated. Black circles 
are the position of states when two nuUclines tangent to each 
other. 
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FIG. 13; Variation of two tangent points uti 2 with respect 

to r. 

1. Uc-i is the "distal" side of the critical cell c. It 
remains in its higher equilibrium since the wave 
back can not pass it. Thus, we can approximate it 
by finding the biggest intersection of the two nuU- 
clines. In the wave back case, since ut-^ is close 
to the higher equilibrium, AUc-i is nearly zero. 
Thus, we determine Wc-i to be 0.9. In the wave 
front case, however, UT2 is small. If we consider 
the minus AUc-i, we determine Wc-i to be 0.8. 

2. Uc+i is the "proximal" side of the critical cell c. It 
switches off before cell c when a wave back comes 
and waits to be excited again by cell c, so at the 
critical time, Uc+i approaches 0. 

3. Figure [T51 shows uti 2 according to Eq. Ob- 
viously, and change little, and can be re- 
garded as the constants 0.7 and 0.17, respectively. 

4. Note that the above approximations are not valid 
when the coupling strength is too large, or ex- 
citability is too weak. 

According to the above approximations, by substitut- 
ing 

MTi = 0.7, Uc-i = 0.9, Uc+i = 
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FIG. 14: The diagram in the parameter plan (Fc-D) repre- 
senting the conditions of collective oscillation. A wave back 
passes the cell, and a wave front reflects, when the parameters 
are above the WB and WF line, respectively. 

into Eq. p^ . and 

UT2 = 0.17, Uc-i ~ 0.8, Uc+i = 

into Eq. (jl6p . we obtain two rough conditions 

b > O.ieSTc ~ 0.7 WB passes, (17) 

D > 0.0398rc + 0.1848 WF generates. (18) 

Clearly, the critical coupling strength increases linearly 
with respect to the active factor F. We draw two lines in 
Fig. [Ml where the WB and WF lines are obtained from 
Eq. ini) and Eq. HH]), respectively. 

In Fig. [lH labels C3, C4 and C5 on the top horizontal 
axis indicate the value of F defined by Eq. ^ for cells 
3, 4 and 5, respectively. Lines WF and WB cross each 
other between C4 and C5. This kind of topology makes 
the A*^ and S*'' cells behave completely different. 

For cell 5, WF is above WB. If the coupling strength 
is between WF and WB, a wave back coming from the 
center can switch the 4*'* cell OFF, but the cell can not 
be switched ON to fire a wave front. This is exactly the 
situation shown in Fig. [5l in which no oscillation takes 
place. On the other hand, for cell 4, WF is below WB. 
Clearly, if the coupling strength allows the wave back to 
suppress the 4*'* cell, the cell will be excited again and 
lead to a wave front. 

The conditions for these two situations depend on 
many other factors, such as the initial conditions, propa- 
gation velocity, cell excitability, and so on. The situation 
is much more complicated than the approximated case 
we have discussed here. Qualitatively, we can conclude 
that the intersection of these two condition lines is the 
origin of self-sustained collective oscillation. 
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V. CONCLUDING REMARKS 

In this paper, we have proposed a scenario for self- 
organized and sclf-sustaincd oscillation in multi-cellular 
biological tissue. In contrast to the usual framework 
based on an oscillatory genetic network, the present sys- 
tem docs not include any sclf-oscillating cells. However, 
by mutual coupling, we can observe collective oscillation 
inside a group of cells, i.e., tissue. Moreover, oscilla- 
tion can manifests itself in several ways, corresponding 
to different coupling strengths. Anti-phase and in-phase 
oscillations at the two boundaries lead to changes in the 
position of oscillation and the oscillating cell population, 
respectively. The birth and death of oscillation resulting 
from variation in the coupling strength were also dis- 
cussed. We also provide a general idea of how the size of 
the cell and population affects the oscillatory behavior. 



Finally, a detailed investigation of the dynamical move- 
ment of the nuUcline provided insight into the mechanism 
of complicated oscillatory phenomena. Although there 
have been several studies on self-oscillatory phenomena 
in spatially discrete systems in the context of mathemat- 
ics and physics, this paper extends these basic ideas to 
spatio-temporal self-organization in a biological system. 
It is of interest to extent our new hypothesis to spatial 
three dimensional systems, i.e., a more realistic model of 
living organism. 

Our observations were based on a numerical simula- 
tion. Future analytical studies inspired by these interest- 
ing phenomenon are needed. At last, but not less impor- 
tant, we are going to cooperate with biologists, in order 
to design corresponding biological experiments and to ex- 
plore more proofs supporting our hypothesis. 
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